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Abstract: Model predictive control (MPC) anticipates future events to take appropriate 
control actions. Nonlinear MPC (NMPC) deals with nonlinear models and/or constraints. 
A Continuation/GMRES Method for NMPC, suggested by T. Ohtsuka in 2004, uses the GMRES 
iterative algorithm to solve a forward difference approximation Ax = 6 of the original NMPG 
equations on every time step. We have previously proposed accelerating the GMRES and 
MINRES convergence by preconditioning the coefficient matrix A. We now suggest simplifying 
the construction of the preconditioner, by approximately solving a forward recursion for the state 
and a backward recursion for the costate, or simply reusing previously computed solutions. 
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1. INTRODUGTION 

Model Predictive Gontrol (MPG) is an optimal control 
technology, which is capable to cope with constrained 
systems and widely used in industry and academia; see, 
e.g., Qin et al. (2003), Camacho et al. (2004), and Griine 
et al. (2011). Nonlinear MPC (NMPC) deals with nonlin¬ 
ear models and/or constraints. Main numerical methods 
applied in NMPC are surveyed by Diehl et al. (2009). 

The Continuation/GMRES by Ohtsuka (2004) is one of 
the real-time numerical methods for NMPC. Ohtsuka’s 
method combines several techniques including replacement 
of inequality constraints by equality constraints, numerical 
elimination of the state, by the forward recursion, and 
the costate, by the backward recursion, and the Krylov 
subspace iterations for solving nonlinear equations via 
parameter continuation. Tanida et al. (2004) have intro¬ 
duced a preconditioned C/GMRES method, however, their 
preconditioner is inefficient. 

Our previous work in Knyazev et al. (2015a) extends 
Ohtsuka’s approach in various ways. The Continuation 
NMPC (CNPMC) method is formulated for a more general 
optimal control model with additional parameters and 
terminal constraints, which allows us solving minimal time 
problems. We also use preconditioners for CNMPC, based 
on an explicit construction of the Jacobian matrices at 
some time steps, improving convergence of the Krylov 
iterations. We propose substituting the MINRES iterative 
solver for GMRES in CNMPC, reducing the memory 
requirements and the arithmetic costs per iteration. 

The present note shows how to reduce the cost of the pre¬ 
conditioning setup, by approximating the Jacobian matrix 
in the Newton iterations. The idea of such an approxi¬ 
mation relies on the observation that most entries of the 


Jacobian weakly depend on small perturbations of the 
state and costate. Most columns of the Jacobian can be 
built from a single instance of the state and costate vari¬ 
ables computed, e.g., during generation of the right-hand 
side of the system solved by the Newton method. Only 
a small number of columns of the Jacobian, specifically, 
responsible for treating the terminal constraints and the 
parameter, is sensitive to changes of the state and costate. 
We recalculate the state and costate corresponding just to 
these sensitive columns. Moreover, for the purpose of the 
preconditioner setup, we can, in addition, compute the 
state and costate on a coarser grid on the horizon with 
subsequent linear interpolation of them at the intermedi¬ 
ate points. We can also use other general techniques for 
fast preconditioner setup, e.g., computation of the state 
and costate variables, as well as the preconditioner and its 
factorization, in a reduced computer precision. Our numer¬ 
ical results demonstrate that the preconditioned GMRES 
and MINRES, where the preconditioner is constructed us¬ 
ing the approximate state and costate variables, converge 
faster, compared to their analogs without preconditioning. 
The paper discusses basic principles of preconditioning, 
and detailed algorithms of computation of the precondi¬ 
tioning schemes are to be reported in our extended paper. 

The rest of the note is as follows. In Section 2, we derive the 
nonlinear equations, which are solved by the continuation 
Newton-Krylov method. Section 3 describes how GMRES 
or MINRES iterations are applied to numerical solution of 
these nonlinear equations. Section 4 presents our main con¬ 
tribution by giving details of the preconditioner construc¬ 
tion, which is based on reusing the previously computed 
and approximated state and costate variables. Section 5 
defines a representative test example; and Section 6 gives 
numerical results illustrating the quality of the method 
with the suggested preconditioner. 




2. DERIVATION OF THE OPTIMALITY 
CONDITIONS 

The MPC approach is based on the prediction by means of 
a finite horizon optimal control problem along a fictitious 
time T € [t,t + T\. Our model finite horizon problem 
consists in choosing the control u{t) and parameter vector 
p, which minimize the performance index J as follows: 

min J, 

u,p 

where 

/•i+T 

J = (l){x{t + T),p) + J L(r,x(r),M(T),p)dr 

subject to the equation for the state dynamics 

^ =/(T,x(r),M(r),p), (I) 

and the equality constraints for the state x and control u 
C{t,x{t),u{t),p) =0, (2) 

V'(a;(t + T),p) = 0. (3) 

The initial value x{T)\r=t for (I) is the state vector x{t) 
of the dynamic system. The control vector u = u(t), 
solving the problem over the prediction horizon, is used 
afterwards as an input to control the system at time t. The 
components of the vector p{t) are parameters of the system 
and do not depend on r. In our minimum-time example 
in Section 5, the scalar parameter p{t) denotes the time to 
destination, and the horizon length is T = p{t). 

The prediction problem stated above is discretized on a 
uniform, for simplicity of presentation, time grid over the 
horizon [t,t -|- T] partitioned into N time steps of size 
At, and the time-continuous vector functions x{t) and 
m(t) are replaced by their sampled values Xi and Ui at 
the grid points n, i = 0,1,... ,N. The integral of the 
performance cost J over the horizon is approximated by 
the rectangular quadrature rule. Equation (I) is integrated 
by the the explicit Euler scheme, which is the simplest 
possible method. We note that more sofisticated one-step 
adaptive schemes can be used as well. The discretized 
optimal control problem is formulated as follows: 

min (j){xN,p) + 

Ui ,p 

subject to 

Xi+i = Xi +f{Ti,Xi,Ui,p)AT, i = 0,l,...,N-l, (4) 

C{Ti,Xi,Ui,p) = 0, i = 0,l,...,N -1, (5) 

'ip{xN,p)=0. (6) 

The necessary optimality conditions for the discretized 
finite horizon problem are the stationarity conditions for 
the discrete Lagrangian function 

N-l 

C{X, U) = (j){xN^p) + ^ L{ri,Xi,Ui,p)Ar 

N-l 

+ Aj[x(t) -a:o] + ^ -Xi+I + f{Ti,Xi,Ui,p)AT\ 

N-l 

+ ^ pJC{Ti,Xi,Ui,p)AT + V^lp{xM,p), 



where X = [xi XN, i = 0,1,... ,N, and U = [ui pi v p]^, 
i = 0,1,... ,N — 1. Here, A is the costate vector, p 
is the Lagrange multiplier vector associated with the 
constraint (5). The terminal constraint (6) is relaxed by 
the aid of the Lagrange multiplier v. 

The necessary optimality conditions are the system of 
nonlinear equations C\^ = 0, Cxi = 0, i = 0,1,... ,N, 
Cu, = 0, = 0, i = 0,1,..., V - 1, C, = 0, Cp, = 0. 

For further convenience, we introduce the Hamiltonian 
function H{t, x, A, u, pL,p) = L{t, x, u,p) -\- f{t, x, u,p) -\- 
pXC{t, X, u,p). 

The optimality conditions are reformulated in terms of 
a mapping F\U,x,t], where the vector U combines the 
control input u, the Lagrange multiplier p, the Lagrange 
multiplier v, and the parameter p, all in one vector: 

U{t) = [ul,... pI ,..., ,p^f. 

The vector argument x in F\U,x,t] denotes the state 
vector at time t, which serves as the initial vector xq in 
the following procedure. 


(1) Starting from the current measured or estimated state 
xq, compute all Xi, i = 0,1..., N, by the forward 
recursion 

Xi+i =Xi + f{Ti,Xi,Ui,p)AT. 

Then starting from 
dd)^ 

Xn = -^{XN,P) + -^{XN,P)V 

compute all costates Xi, i = N,...,1,0, by the 
backward recursion 

Ai — Aj-|_i T ^ (tj, , Ai-|-i, 

(2) Using just obtained Xi and Xi, calculate the vector 

r ou^ 

— — {to,xo,Xi,uo,po,p)At 


F[U,x,t] = 


-{Ti,Xi,X^+i,Ui,Pi,p)AT 


-(tat-i, xat-i, Aat, un-1, Pn-i,p)At 
C{to,xo,uo,p)At 


[Ti.Xi. Uj.pjlAT 


C{tn-i,xn-i,un-i,p)At 

'tp{xN,P) 

ddX dilX 

-^i^N,p) + -^{xN,p)y 

-^-1 qpjT 

+ T! “o— {n,Xi,X,+i,Ui,pi,p)AT 


The equation with respect to the unknown vector U(t) 

F[Uit),xit),t]=0 (7) 

gives the required necessary optimality conditions. 



3. NUMERICAL ALGORITHM 


4. PRECONDITIONING 


The controlled system is sampled on a uniform time grid 
tj = jAt, j = 0,1,.... Solution of equation (7) must 
be found at each time step tj in real time, which is a 
challenging part of implementation of NMPC. 

Let us denote xj = x{tj), Uj = U{tj), and rewrite the 
equation F\Uj,Xj,tj] = 0 equivalently in the form 

F\Uj , Xj , t] — F\Uj-i^ Xj,tj] = bj, 

where 

bj = —F[Uj-i,Xj,tj]. (8) 

Using a small h, which may be different from At and At, 
we introduce the operator 

aj{V) = {F[Uj-i + hV,Xj,tj] - F[Uj-i,Xj,tj])/h. (9) 

We note that the equation F[Uj,Xj,tj] = 0 is equivalent to 
the equation aj{AUj/h) = bj/h, where AUj = Uj — Uj-i. 

Let us denote the fc-th column of the mxm identity matrix 
by Cfc, where m is the dimension of the vector U, and define 
an m X m matrix Aj with the columns AjCk = aj{ek), 
k = 1, ... ,m. The matrix Aj is an 0{h) approximation of 
the Jacobian matrix Flf[Uj-l,Xj,tj]^ which is symmetric. 

Suppose that an approximate solution Uq to the initial 
equation F[Uo, xq, to] = 0 is available. The first block entry 
of Uq is taken as the control uq at the state xq. The next 
state xi = x{ti) is either sensor estimated or computed by 
the formula xi = xq + f{to,xo,uo)At; cf. (I). 

At the time tj, j > 1, we have the state Xj and the vector 
Uj-i from the previous time tj-i. Our goal is to solve the 
following equation with respect to V: 

ajiV) = b,fh. (10) 

Then we can set Uj = Uj-i + hV and choose the first block 
component of Uj as the control Uj . The next system state 
Xj+i = x{tj+i) is either sensor estimated or computed by 
the formula Xj+i = Xj + f{tj,Xj,Uj)At. 

A direct way to solve (10) is generating the matrix Aj and 
then solving the system of linear equations AjAUj = bj; 
e.g., by the Gaussian elimination. 

A less expensive alternative is solving (10) by the GMRES 
method, where the operator aj{V) is used without explicit 
construction of the matrix Aj (cf., Kelly (1995); Ohtsuka 
(2004)). Some results on convergence of GMRES in the 
nonlinear case can be found in Brown et al. (2008). 

We recall that, for a given system of linear equations 
Ax = b and initial approximation xg, GMRES con¬ 
structs orthonormal bases of the Krylov subspaces /C„ = 
spanjro, Aro,..., A”“^ro}, n = 1 , 2 ,..., given by the 
columns of matrices such that AQn = Qn+iHn with 
the upper Hessenberg matrices and then searches for 
approximations to the solution x in the form Xn = QuUn, 
where ?/„ = argmin|| - ^Ib- 

A more efficient variant of GMRES, called MINRES, may 
be applied when the matrix A is symmetric, and the 
preconditioner is symmetric positive definite. Using the 
MINRES iteration in Ohtsuka’s approach is mentioned in 
Knyazev et al. (2015a). 


The convergence of GMRES can be accelerated by pre¬ 
conditioning. A matrix M that is close to the matrix A 
and such that computing M~^r for an arbitrary vector 
r is relatively easy, is referred to as a preconditioner. 
The preconditioning for the system of linear equations 
Ax = b with the preconditioner M formally replaces the 
original system Ax = b with the equivalent precondi¬ 
tioned linear system M~^Ax = M~^b. If the condition 
number ||M“^A|| ||of the matrix M~^A is small, 
convergence of iterative Krylov-based solvers for the pre¬ 
conditioned system can be fast. However, in general, the 
convergence speed of, e.g., the preconditioned GMRES is 
not necessarily determined by the condition number alone. 

A typical implementation of the preconditioned GMRES 
is given below. The unpreconditioned GMRES is the same 
algorithm but with M = I, where I is the identity matrix. 
We denote by the submatrix of H with the 

entries H^j such that ii < i < i 2 and ji < j < j 2 - 

Algorithm Preconditioned GMRES(fcmax) 

Input: a{v), b, xg, fc^ax, M 
Output: Solution x of a(x) = b 

r = b- a{xg), z = M-'^r, (3 = \\z\\ 2 , Wi = z/f3 
for A: = 1,..., fcmax do 
r = a{vk), z = M~^r 
Hl:k,k = [vi, . . . ,Wfc ]^2 
z = z - [ui, .. .,Vk]Hi,k,k 
Hk+i,k = ||2:||2 
Vk+l = z/\\z \\2 
end for 

y = arg minj^||Hi:fc„„+i,i:fc„„ 2 / - [/3, 0,... ,0]^||2 
x = xo + [vi,...,Vk^,Jy 

In Knyazev et al. (2015a), the matrix Aj is exactly com¬ 
puted at some time instances tj and used as a precondi¬ 
tioner in a number of subsequent time instances tj, tj+i, 
..., tj+jp. In the present note, we propose to use a close 
approximation to Aj, which needs much less arithmetic 
operations for its setup. Construction of such approxima¬ 
tions Mj is the main result of this note. 

We recall that computation of the fc-th column of Aj 
requires computation of all states x(Ti) and costates X{Ti) 
for the parameters stored in the vector Uj-i + hck- Is it 
possible to replace them by x{Ti) and A(ri) computed for 
the parameters stored in the vector Uj-il The answer is 
yes, for the indices k = 1,... ,m — I, where I is the sum of 
dimensions of ip and p. These k indices correspond to the 
terms containing the factor At in the Lagrangian C. 

The first m — l columns (and rows, since the preconditioner 
Mj is symmetric) are calculated by the same formulas as 
those in Aj, but with the values x(Ti) and A(Ti) computed 
only once for the parameters stored in the vector Uj-i, 
i.e., when computing the vector bj. Thus, the setup of Mj 
computes the states x{Ti) and costates \{Ti) only I times 
instead of m times as for the matrix Aj. It is this reduction 
of computing time that makes the preconditioner Mj more 
efficient, especially in cases where dimension of the state 
space is very large. 



The preconditioner Mj is obtained from Aj by neglect¬ 
ing the derivatives dxi^/dui^, dXi^/dui^, dxi^/d^i^ and 
dXi^/dfii^. Therefore, the difference Aj — Mj is of order 
0 (At) since dxi^/dui^ = 0(Ar), dXi^/dui^ = 0(Ar), 
dxi^/djii^ = 0 and dXi^/d^i^ = 0 {At). 

The preconditioner application M~^r requires the LU 
factorization M = LU, which is computed by the Gaussian 
elimination. Then the vector M~^r = U~^{L~^r) is ob¬ 
tained by performing back-substitutions for the triangular 
factors L and U. Further acceleration of the preconditioner 
setup is possible by faster computation of the LU factoriza¬ 
tion. For example, when computation with lower number 
of bits is cheaper than computation with the standard 
precision, the preconditioner Mj and its LU factorization 
may be computed in lower precision. 

Another way of reduction of the arithmetical work in the 
preconditioner setup is the computation of the states x 
and costates A with the double step 2 At thus halving the 
arithmetical cost and memory storage. The intermediate 
values of x and A are then obtained from the computed 
values by simple linear interpolation. 

5. EXAMPLE 


We consider a test nonlinear problem, which describes 
the minimum-time motion from a state {xo,yo) to a state 
{xf,yf) with an inequality constrained control: 


State vector x = 


and input control u = 


u 

Ud 


Parameter variable p = tf — t, where tf denotes the 
arrival time at the terminal state {xf,yf). 

Nonlinear dynamics is governed by the system of 
ordinary differential equations 

(Ax ■ 


X = f{x,u,p) = 


{Ax ■ 


B) cosu 
■ B) sinw 


Constraint: C{x,u,p) = [u — Cuf + u\ — r'^, 


= 0 , 


where Cu = cq + ci sin(a;t) and Ud is a slack variable, 
i.e., the control u always stays within the sinusoidal 
band c„ - r„ < u < c„ -I- r„). 

"x-xf 


Terminal constraints: 'ip{x,p) = 


y-vf 


= 0 (the 


state should pass through the point {xf ,yf) a,t t = tf) 
Objective function on the horizon interval [t,tf]' 

pt+p 

J = (j){x,p)+ / L{x,u,p)dt, 

where 


4>{x,p) =p, L{x,u,p) = -WdUd 

(the state should arrive at {xf,yf) in the shortest 
time; the function L serves to stabilize the slack 
variable Ud) 

• Constants: A = B = 1, xq = yo = 0, Xf = yf = 1, 
Co = 0.8, Cl = 0.3, uj = 10, r„ = 0.2, Wd = 0.005. 


The horizon interval [t,tf] is parametrized by the affine 
mapping t ^ t + rp with t G [0,1]. 

The components of the corresponding discretized problem 
on the horizon are given below: 

• At = 1/A, Tj = iAr, Cm = co + ci sm{uj{t + np)); 


the participating variables are the state 


, the 


costate 


Ai,i 

X2,i 


, the control 


^di 


, the Lagrange 


multipliers pi and 


Vl 

V2 


, the parameter p; 


the state is governed by the model equation 


f x^+l = x^ + At [p {Axi + B) cosm*] , 
\ j/j+i =yiL At [p {Ax^ -I- B) sin ui ], 


where f = 0,1,..., A — 1; 

• the costate is determined by the backward recursion 
(Ai,Ar = Vl, A2,Ar = V 2 ) 


Ai,i 

A2,i 


— Ai_i+i 

-I- At [pA{cosuiXij+i 

= A2,i+1, 


-|-sinMA2.*-i-i)], 


where f = A — 1,A — 2, ...,0; 

• the equation F{U,xo,yo,t) = 0, where 


U = [uo, Ud,0, ■ • ■ , UN-l,Ud,N-l, 

Po,---,Pn-i,vi,V2,p], 

has the following rows from the top to bottom: 

At [p{Axt + B) (-sinuiAi,i-i-i + cosmA 2 ,*-i-i) 

2 {tLi Cui) Pi\ — 0 

At [2piUdi - Wdp] = 0 
{ At [{ui - Cui)'^ + uli -rl] =0 

J XJV — Xr = 0 
\yN -yr = 0 
Af-l 

At[^ {Axt -I- B){cosuiXij+i + sinMA2,*-i-i) 

i—0 

- 2{ui - Cui)piCi cos{ui{t + np))ujTi 

-WdUdi] -1-1 = 0. 

Let us compare the computation costs of the matrices Aj 
and Mj for this example. We do not take into account the 
computation of the right-hand side bj = —F[Uj-i,Xj,tj] 
because it is a necessary cost. Computation of the matrix 
Aj requires 3A -|- 3 evaluations of the vector F[Uj-i + 
hV,Xj,tj], where A is the number of grid points on the 
prediction horizon. Setup of Mj requires only 3 evaluations 
of F\Uj-i + hV,Xj,tj], which is A -|- 1 times faster. 


6. NUMERICAL RESULTS 


In our numerical experiments, the weakly nonlinear sys¬ 
tem (10) for the test problem from Section 5 is solved 
by the GMRES and MINRES iterations. The number of 
evaluations of the vector a{V) at each time tj does not 
exceed an a priori chosen constant /cmax = 20. In other 
words, the maximum number of GMRES or MINRES 
iterations is less or equal /cmax- The error tolerance in 
GMRES and MINRES is tol = I0“®. The number of grid 
points on the horizon is A = 50, the sampling time of 
simulation is At = 1/500, and h = 10“®. 

The preconditioners are set up at the time instances Up, 
where tp = 0.2 is the period, and I = 0,1,.... After each 
setup, the same preconditioner is applied until next setup. 
Preconditioners for MINRES must be symmetric positive 
definite and are built here as the absolute value of Mj, i.e., 



if Mj = UY.V'^ is the singular value decomposition, then 
\Mj\ = see Vecharinski et al. (2013). 

Figure 1 shows the computed trajectory for the test 
example. Figure 2 shows the optimal control by the MFC 
approach using the preconditioned GMRES. Figure 3 
displays ||E ||2 and the GMRES residuals. 

The number of iterations of preconditioned GMRES is 
displayed in Eigure 4. Eor comparison, we show the number 
of iterations of the preconditioned MINRES in Figure 5. 

Figure 6 displays ||E ||2 and the 2-norm of the residual 
after iterations of GMRES without preconditioning. The 
corresponding number of iterations of GMRES without 
preconditioning is shown in Eigure 7. The number of 
iterations of MINRES without preconditioning is shown 
in Figure 8. 

Effect of preconditioning is seen when comparing Eigures 4 
and 7 for GMRES and Figures 5 and 8 for MINRES. 
Preconditioning of GMRES reduces the number of iter¬ 
ations by factor 1.2. Preconditioning of MINRES reduces 
the number of iterations by factor 1.4. 

The number of iterations does not necessarily account for 
the additional complexity that preconditioning brings to 
the on-line algorithm. However, the computation time is 
machine and implementation dependent, while our tests 
are done in MATLAB on a generic computer. Specific 
implementations on dedicated computer chips for on-line 
controllers is a topic of future work. 

GONCLUSIONS 

We have found a new efficient preconditioner Mj, which 
approximates the Jacobian matrix Aj of the mapping F 
defining equation (7). Gomputation of Mj is 0{N) times 
faster than that of Aj, where N is the number of grid 
points on the prediction horizon. 

The preconditioner Mj can be very efficient for the NMPG 
problems, where dimension of the state space is large, for 
example, in the control of dynamic systems described by 
partial differential equations. 

Other useful techniques for accelerating the preconditioner 
setup include computation of the matrices Mj and their 
LU factorizations in lower precision, computation of the 
state and costate on a coarse grid over the horizon and 
linear interpolation of the computed values on the fine grid. 
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Fig. 1. Trajectory by NPMC using the preconditioned 
GMRES 



Fig. 2. NMPC control u using the preconditioned GMRES 










Fig. 3. Preconditioned GMRES, fcmax = 20 
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Fig. 6. GMRES without preconditioning, fcmax = 20 



Fig. 4. Preconditioned GMRES, fcmax = 20 
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Eig. 7. GMRES without preconditioning, fcmax = 20 




Fig. 5. Preconditioned MINRES, fcmax = 20 


Eig. 8. MINRES without preconditioning, fcmax = 20 












































































